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(57) Abstract: A method of image velocity estimation in image processing which uses a block matching technique in which a 
similarity measure is used to calculate the similarity between blocks in successive images. The similarity measure is used to calculate 
a probability density function of candidate velocities. The calculation is on the basis of an exponential function of the similarity in 
which the similarity is multiplied by a parameter whose value is independent of position in the frame. The candidate velocities 
are thresholded to exclude those having a low probability. The value of the parameter and threshold are optimised together by 
coregistering all frames to the first frame, calculating the registration error, and varying them to minimise the registration error. The 
similarity measure is normalised with respect to the size of the block, for example by dividing it by the number of image samples 
in the blocks being compared. The similarity measure used may be the CD2-bii similarity measure in which the mean and standard 
deviation of the two blocks being compared are adjusted to be the same before calculation of the similarity. This makes the similarity 
measure particularly suitable for ultrasound images. Further, block matching may be conducted across three frames of the sequence 
by comparing the intensities in blocks in the first and third, and second and third of the frames and finding the block in the third 
frame which best matches ihe block in the second frame and that block's corresponding position in the first frame. 
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The present invention relates to image processing, and in particular to 
improving the estimation of image velocity in a series of image frames. 

There are many imaging situations in which a subject in an image is in motion 
and it is desired to track or measure the movement of the subject from frame to 
frame. This movement is known as optical flow or image velocity. Such estimation 
or measurement of image velocity may be done, for example, to improve the 
efficiency of encoding the image, or to allow enhancement of the display of, or 
measurement of, the movement of some particular tracked part of the image to assist 
an observer trying to interpret the image. Many techniques have been proposed and 
used for image velocity estimation and one of the basic techniques is known as block 
matching. In block matching, blocks of pixels are defined in a first frame and the 
aim is then to identify the position of those blocks in a second subsequent.frame. 
One approach is to compare the intensities of the pixels in the block in the first frame 
with successive, displaced candidate blocks in the second frame using a similarity 
measure, such as the sum of square differences. The block in the second frame 
which gives the minimum of the sum of square differences (or gives the best match 
with whichever similarity measure is chosen) is taken to be the same block displaced 
by movement of the subject. Repeating the process for successive blocks in the first 
image frame gives an estimate for the subject motion at each position in the image 
(the image velocity field). 

Figure 1 schematically illustrates the idea. Two frames are shown, frame 1 
and frame 2. These may be, but are not necessarily, successive frames in a sequence. 
Frame 1 is divided up into square blocks of pixels having a side length of (2 n + 1) 
pixels, ie. from -n to +n about a central pixel (x, y) in each block. One block W c is 
illustrated in Fig. 1 . A search window W s is defined in the second frame around the 
position of the corresponding central pixel (x, y) in the second frame. As illustrated 
in Fig. 1 it is a square search region of side length (2 N + 1) pixels. The intensities of 
the block W c of pixels in frame 1 are then compared at all possible positions of the 



WO 2004/052016 PCT/GB2003/005047 



-2- 

block in the search window W s . So, for example, the first comparison is made with 
the corresponding (2 n + 1) by (2 n + 1) block in the top left hand corner of the search 
window W s , and then with such a block displaced one pixel to the right, and then a 
block displaced two pixels to the right and so on until the end of the search window 
5 is reached. The procedure is then repeated for a row of candidate blocks displaced 
one pixel down in the search window from the first row, and so on until the bottom 
of the search window is reached. The similarity measure may, for example, be a sum 
of square differences:- 

n n 

10 £ c («,v)=n |/(* +**y+ J, t)-I(x+u+i,y+v+j,t + l)] 2 (1) 

i=~nj--n 

for each value of (w , v) for -N^^v^and where i and j index through the block W c 
centered on the pixel (x , y) in the x and y directions respectively, and u and v are the 
different values of displacement which index over the search window W 5 . The 
values u and v can, given the time difference between the frames, be regarded as a 

1 5 velocity. This gives a value of E c for each estimated displacement. The estimated 
displacement with the minimum E c is often taken as the actual displacement of the 
block. This is repeated for all positions in frame 1 to give a velocity field, and then 
for all frames in the sequence. Different similarity measures may, tff course, be used. 
Also, it is not always necessary to conduct the block matching on all frames of the 

20 sequence, or for all pixels or blocks in each frame. The block W c may subsample the 
pixels in the frame and the candidate displacements u and v may be indexed by more 
than one pixel. Thus the searching may be at different resolutions and scales. 
Sometimes a multi-scale and/or multi-resolution approach may be used in which 
block matching is first performed at a coarse resolution or large scale, and 

25 subsequently at successively finer resolutions, using the previously calculated 
velocity values to reduce the amount of searching required at finer resolutions. 

Medical images present many difficulties in image processing because of the 
typically high level of noise found in them. For example, the tracking of cardiac 
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waUs in cardiac ultrasound images is difficult because of the high level of noise 
found in ultrasound images and also because of the nature of the cardiac motion. In 
particular, the cardiac motion varies during the cardiac cycle. Various ways of 
identifying and tracking cardiac walls in echocardiograms have been proposed in WO 
5 01/16886 and WO 02/43004, but it is a difficult task in which there is room for 
improvement. 

A development of the block matching technique as described above has been 
proposed by A. Singh, "Image-flow computation: An estimation-theoretic framework 
and a unified perspective," CVGIP: Image understanding, vol. 65, no. 2, pp. 152-177, 

10 1 992 which is incorporated herein by reference. In this approach both conservation 
information, e.g. from a block matching technique as described above, and 
neighbourhood information (i.e. looking at the velocities of surrounding pixels) are 
combined with weights based on estimates of their associated errors. Thus in a first 
step based on conservation information the similarity values E c are used in a 

15 probability mass function to calculate a response R c whose value at each position in 
the search window represents the likelihood of the corresponding displacement. The 
probability mass function is given by 

1 

R c (u,v)=-exp(-kE c (u t v)) 9 (2) 
where Z is defined such that all probabilities sum to unity i.e.:- 

20 

EN ^ , . A 

In the function for the response the parameter k is chosen at each position 
such that the maximum response in the search window is close to unity (0.95 before 
normalisation) for computational reasons. The expected value of the velocity is then 
25 found by multiplying each candidate value by its probability and summing the 
results:- 
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U, 



N 



N 



uR c (u,v) (4) 



cc 




N 



vR c (u,v) 



(5) 



10 



15 



Another velocity estimate may be obtained by the use of neighbourhood 
information. In other words, the velocity at each pixel is unlikely to be completely 
independent of the velocity of its neighbours. Thus, assuming that the velocity of 
each pixel in a small neighbourhood W p has been estimated, the velocity estimates 
for each pixel can be refined by using the velocity of its neighbouring pixels. Clearly 
it is more likely that the velocities of closer neighbours are more relevant than pixels 
which are further away. Therefore weights are assigned to velocities calculated for 
the neighbouring pixels, and the weights drop with increasing distance from the 
central pixel (a 2-D Gaussian mask in the window W p of size (2w+l)(2w+l) is used). 
These weights can be interpreted as a probability mass function R n = (u i3 v,) 
where X i€W P Rn ( Ui > v< ) = 1 ( i is an index for pixels in W p ) in a uv space. 
Now, the velocity estimate U = ( u, v ) for the central pixel using neighbourhood 
information can be calculated as: 



An important aspect of the Singh approach is that a covariance matrix is 
calculated for each velocity estimate, for both the estimates based on conservation 
information and the estimates based on neighbourhood information. These 
covariance matrices can be used to calculate errors which are used as weights when 
combining the two estimates together to give a fused, optimal estimate. 




) (6) 




(7) 
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The covariance matrix corresponding to the estimate U a is given by:- 



£l_„Z u a ) 2 KM El.„El.v,("- u cc)(y- vjRX^v) 



(8) 



The covariance matrix corresponding to the neighbourhood estimate U is as 



5 follows:- 




E, <6 r,fa " "Xv, - vK( M ,,v,) " v) 2 ^(«„ v,-) 



(9) 



Thus these steps give two estimates of velocity, U cc and £/ , from 
conservation and neighbourhood information respectively, each with a covariance 
1 0 matrix representing their error. An estimate U of velocity that takes both 

conservation information and neighbourhood information into account can now be 

computed. The distance of this new estimate from U , weighted appropriately by the 
corresponding covariance matrix, represents the error in satisfying neighbourhood 
information. This can be termed neighbourhood error. Similarly the distance of this 
15 estimate from weighted by its covariance matrix, represents the error in 

satisfying conservation information. This may be termed conservation error. The 
sum of neighbourhood and conservation errors represents the squared error of the 
fused velocity estimate:- 



* 2 = (u- u„) r sg{u- uj+ (u- u) T s;\u- u) (io) 



20 



The optimal value of velocity is that value which minimises this error and can 
be obtained by setting the gradient of the error with respect to U equal to zero 



giving:- 
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= + s l Y l [sgu„ + s?u\ (ii) 

Because the values of U and iS^ are derived on the assumption that the 

velocity of each pixel of the neighbourhood is known in advance, in practice equation 
(1 1) is solved in an iterative process (via Gauss-Seidel relaxation) with the initial 
5 values of the velocity at each pixel being taken from the conservation information 
alone. Thus:- 

\k 1 = [s;^s^[s;;u cc+ s;^} (12) 

where the superscript m refers to the iteration number. Iteration continues 
until the difference between two successive values of U op is smaller than a specified 
10 value. 

While this technique usefully combines conservation and neighbourhood 
information, and does so in a probabilistic way, it does not always work well in 
practice, particularly with noisy images of the type found in medical imaging and 
ultrasound imaging in general. 

1 5 Another common problem in image velocity estimation using matching 

techniques is known as the multi-modal response (i.e. due to the well-known aperture 
problem, for example, or mismatching especially when the size of the search window 
is large). A common way to reduce the effect of the multi-modal response is to 
compare the intensities in three frames, rather than just two as explained above. So 

20 the similarity between blocks W c in two frames x { and y { is found, and between two 
blocks W c my i and z u as shown in Figure 2 of the drawings. In the two frame 
comparison the intensities in a block W c in one frame x t at time t are compared with 
the intensities in a corresponding block displaced by the candidate velocity (w,v) in 
the next frame y ( at time t+1 for all values of (u, v) in the search window W s . In the 

25 three frame approach, the intensities in the block W c are also compared with the 

intensities in the block displaced by (2u, 2v) in the next-but-one frame z x at time t+2, 
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again for values of (w, v) in the search window W s . In the case of using sum-of- 
square differences as the similarity measure this can be written as: 



£ c («,v)=II 



i=-nj=-n 



[l(x + i 9 y + j, t)-I( x +u+i,y+v+j,t+ l)f 

+ (13) 

[l(x + Uy + y\0 - i(jc + 2u + /, y + 2v + y, f + 2)] 2 



5 where the first term is comparing blocks in the frames at t and t+1 separated 

by a displacement (u, v) and the second term is comparing blocks in the frames at t 
and t+2 separated by twice that, i.e. (2u, 2v). This amounts to assuming that the 
velocity is constant across three frames of the sequence. In other words, for three 
frames at times t, t+1 and t+2, it is assumed that the displacements between t and t+1 

10 are the same as the displacements between t+1 and t+2. This assumption is 

reasonable for high frame rate sequences, but is poor for low frame rate sequences, 
such as are encountered in some medical imaging techniques, including some 
ultrasound imaging modalities. 

The present invention is concerned with improvements to block matching 

1 5 which are particularly effective for medical images, especially ultrasound images, 
which are inherently noisy. 

A first aspect of the invention provides a method of processing a sequence of 
image frames to estimate image velocity through the sequence comprising: 

block matching using a similarity measure by comparing the intensities in 

20 image blocks in two frames of the sequence and calculating the similarity between 
the said blocks on the basis of their intensities, calculating from the similarity a 
probability measure that the two compared blocks are the same, and estimating the 
image velocity based on the probability measure, wherein the probability measure is 
calculated using a parametric function of the similarity which is independent of 

25 position in the image frames. 
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Preferably the parameters of the parametric function are independent of 
position in the image frames. The function may be a monotonic, e.g. exponential, 
function of the similarity, in which the similarity is multiplied by a positionally 
invariant parameter. The parameters may be optimised by coregistering the frames in 
5 the sequence on the basis of the calculated image velocity, calculating a registration 
error and varying at least one of the parameters to minimise the registration error. 
The registration error may be calculated from the difference of the intensities in the 
coregistered frames, for example the sum of the squares of the differences. 

Thus in the particular example of the approach proposed by Singh the value 

10 of parameter k is set for each position (so that the maximum response in the search 
window is close to unity), meaning that k varies from position to position over the 
frame. However, with this aspect of the present invention the value of k is fixed over 
the frame - it does not vary from position to position within the frame. It should be 
noted that because k is used in a highly non-linear (exponential) function in 

15 calculating the response (probability), the velocity and error estimates are not 

uniform, because variations in the value of k have a large effect. With this aspect of 
the present invention, on the other hand, k is constant for all pixels in the image, so 
the processing is uniform across the image and from frame to frame. 

The value of k may be optimised, as mentioned, for example by registering all 

20 frames in the sequence to the first frame, i.e. using the calculated image velocity to 
adjust the image position to cancel the motion - which if the motion correction were 
perfect would result in the images in each frame registering perfectly, and calculating 
the registration error - e.g. by calculating the sum of square differences of the 
intensities. The value of is chosen which gives the minimum registration error. 

25 The calculated similarity may be normalised by dividing it by the number of 

pixels in the block, or the number of image samples used in the block (if the image is 
being sub-sampled). 

Thus, again in the particular example above, the value of k in equation (2) 
above for R c may be replaced by k/(2n+lf. This means that the value of k does not 

30 need to be changed if the block size is changed. In particular, it does not need to be 
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re-optimised, so that once it has been optimised for a given application (e.g. breast 
ultrasound) using one frame sequence at one scale and resolution, the same value of k 
may be used for the same application on other sequences at other scales and 
resolutions. 

5 The probability measure may be thresholded such that motions in the image 

velocity having a probability less than a certain threshold are ignored. The threshold 
may be optimised by the same process as used for optimisation of the parameter k 
above, i.e. by coregistering the frames in the sequence on the basis of the calculated 
image velocity, calculating a registration error and varying the threshold to minimise 

1 0 registration error. The threshold may be positionally independent. 

A second aspect of the invention relates to the similarity measure used in 
image velocity estimation and provides that the intensities in the blocks W c in the 
frames being compared are normalised to have the same mean and standard 
deviations before the similarity is calculated. The similarity measure may be the CD 2 

1 5 similarity measure (rather than the sum of square differences of Singh), which is 
particularly suited to ultrasound images (see B. Cohen and I. Dinstein, "New 
maximum likelihood motion estimation schemes for noisy ultrasound images", 
Pattern Recognition 35 (2002), pp 455-463). 

A third aspect of the invention modifies the approach of Singh to avoiding 

20 multi-modal responses by assuming that the observed moving tissue conserves its 
statistical behaviour through time (at least for three to four consecutive frames), 
rather than assuming a constant velocity between three frames. 

This aspect of the invention provides for block matching across three frames 
of the sequence by comparing the intensities in blocks in the first and third and the 

25 second and third of the three frames, and calculating the similarity on the basis of the 
compared intensities. 

The blocks in the first and second frames are preferably blocks calculated as 
corresponding to each other on the basis of a previous image velocity estimate (i.e. 
the image velocity estimate emerging from processing preceding frames). 

30 Thus the method may comprise defining for each block in the second frame a 
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search window encompassing several blocks in the third frame, and calculating the 
similarity of each block in the search window to the said block in the second frame 
and to the corresponding position of that said block in the first frame (as deduced 
from the previous image velocity estimate). Thus this avoids assuming that the 

5 velocity remains the same through the three frames. It is therefore suited to image 
frame sequences having a relatively low frame rate, where the assumption of constant 
velocity does not tend to hold. 

The different aspects of the invention may advantageously be combined 
together, e.g. in an overall scheme similar to that of Singh. Thus, as in the Singh 

0 approach the estimated image velocity using the technique above may be obtained by 
summing over the search window the values of each candidate displacement 
multiplied by the probability measure corresponding to that displacement. Further, 
the estimate may be refined by modifying it using the estimated image velocity of 
surrounding positions - so-called neighbourhood information. 

5 The techniques of the invention are particularly suitable for noisy image 

sequences such as medical images, especially ultrasound images. 

The invention also provides apparatus for processing images in accordance 
with the methods defined above. The invention may be embodied as a computer 
program, for example encoded on a storage medium, which executes the method 

0 when run on a suitably programmed computer. 

The invention will be further described by way of example, with reference to 
the accompanying drawings in which:- 

Fig. 1 illustrates schematically a block matching process; 
5 Fig. 2 illustrates schematically a similarity measure calculation using a 

constant velocity assumption for three frames; 

Fig. 3 illustrates a similarity measure calculation using the assumption of 
statistical conservation of moving tissue for three frames; 

Fig. 4 is a flow diagram of an optimisation process used in one embodiment 
0 of the invention; 
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Fig. 5 illustrates the overall process of one embodiment of the invention; and 
Fig. 6 illustrates the optimisation of k and 7 for a breast ultrasound image 
sequence. 

5 Given a sequence of image frames in which it is desired to calculate the 

image velocity, the first aspect of the invention concerns the similarity measure used, 
i.e. the calculation of E e (u, v). While the image processing algorithm proposed by 
Singh uses the sum of square differences as a similarity measure, other similarity 
measures such as CD 2 and normalised crossed correlation (NCC) are known. In this 
10 embodiment a modified version of the CD 2 similarity measure is used. Using the 
CD 2 similarity measure the most likely value of the velocity is defined as:- 

vT = max E x v ~ y*J - ln(exp(2(% - y tJ ) + 1) (14) 
v, >' 

where here i refers to the block, y indexes the pixels in the block, there are 2n+l 
pixels in the block, and x fj and y tJ are the intensities in the two blocks being 
15 compared. 

This similarity measure is better for ultrasound images than others such as 
sum-of-square differences or normalised cross-correlation because it takes into 
account the fact that the noise in an ultrasound image is multiplicative Rayleigh 
noise, and that displayed ultrasound images are log-compressed. However it assumes 

20 that the noise distribution in both of the blocks W c is the same and this assumption is 
not correct for ultrasound images. The attenuation of the ultrasound waves 
introduces inhomogeneities in the image of homogeneous tissue. The time gain and 
the lateral gain compensations (compensating respectively for the effects that deeper 
tissue appears dimmer and for intensity variations across the beam) which are tissue 

25 independent and generally constant for a given location during the acquisition, do not 
compensate fully for the attenuation. Thus in this embodiment of the invention an 
intensity normalisation is conducted before calculation of the CD 2 similarity measure. 
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This is achieved by making sure that the two blocks W c of data have at least the same 
mean and variance. In more detail, the original intensity values and ja. above are 

transformed into new values of x tj and y fJ by subtracting the mean and dividing by 

the standard deviation (square root of the variance) of the intensity values in the 
5 block. This gives a new similarity measure which can be denoted CD 2 . bis as follows:- 

2n+l 

Ef^ to(exp(2(x, - y 9 )) + 1) (15) 

This is the similarity measure used in this embodiment to calculate the values 
of 2? a used. 

To avoid multi-modal responses, the similarity measure may be calculated 

1 0 over three consecutive frames. However, rather than making the normal constant 
velocity assumption as mentioned above and described in relation to Figure 2, which 
results in the similarity measure being based on comparing the first frame at time t 
with the next frame at time t+1 and the third frame at t+2, instead the result of 
calculating the velocities between the preceding frame at time t-1 and the current 

1 5 frame at time t are used. Given a block in frame x { at time t, which is compared to 
blocks in the search window f^in frames at the time t+1, the position of that block 
in the preceding frame at time t-1 (denoted has already been calculated and so its 
position can be denoted (x-u^ y-v 0 ) in the preceding frame where (u^ v 0 ) was the 
result of the preceding velocity (image velocity) calculation. Thus in the three frame 

20 approach in this embodiment of the invention the intensities of each candidate block 
in the search window W s are compared with the intensities of the block at (x, y) in the 
frame x t at time t, and also with the calculated position (x-Up y-v 0 ) of that block in the 
frame o i at time t-1. A value of E is calculated for each comparison (of and and 
o l andj/j) and the values are summed. 

25 This is illustrated schematically in Figure 3. The approach is applicable 

whatever similarity measure is used to compare the intensities. In the case of the 
sum-of-square differences, the new similarity measure becomes:- 
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[l(x - u 0 + i , y - v Q + ; , t - 1) - I{x + u + i , y + v + ; , f + l)f . 



(16) 



([l(x+ i 9 y+ ;,/)- /(* + v+ j,t + 1)] 2 



where the first term compares intensities in frames and#, i.e. at times t-1 and t+1, 
and the second term compares intensities between frames x i and j/„ i.e. at times t and 
5 t+1. 

In the case of CD^ similarity measure defined above, the calculation of E 
over three frames becomes:- 



2w+l 



£ °y - Jy - ln(exp(2(5; y - y tJ )) + 1) 



(17) 



10 or in more detail: 



7(x- 1^+ v 0 + 1)- 7(x + r/+ /,^+ v + y,f + 1)- 
v ^ ^ [ln(exp(2[7(x- ^ + v 0 + j\t- 1)- 7(x + «+ f f jr+ v+ y.f + 1)])+ 

^ln(exp(2[7(x+ i, j/+ j t t)~ 7(x+ u+ v+ j\t + 1)])+ l) ; 



(18) 



Here / represents the intensity data I transformed as detailed above (but only, of 
course, within the interesting block, not for the whole image). 

This avoids the assumption that the velocity is the same over the three frames. 
1 5 Instead it looks for the best match in frame y i to the block in x t and the calculated 
previous position of the block (in o^. It improves the matching process especially at 
low frame rates, e.g. of 20-30 Hz. This makes it particularly useful in the case of 
contrast echocardiography, abdominal imaging, tissue Doppler and real-time 3D- 
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imaging, where low frame rates are typical. 



Having established the new similarity measure, the next stage is to 
5 calculate a probability mass function from the similarity measure. In the Singh 

approach this was by equation (2) above. As discussed above, that involved setting a 
value of k for each position in the frame such that the maximum response in the 
search window was close to unity. However, in this embodiment of the invention the 
value of k is, instead, set to be the same for all positions in the frame and all frames 
10 in the sequence. The value of k is found in an optimisation approach which will be 
described below. Given the value of k the probability mass function for this 
embodiment is given by 

( i \ 



1 

#A v;= — exp 



U2#i+1> 



2. 



(Ec(u 9 v)-m) 



(19) 



where m is the maximum of the similarity measure in the search window W s (i.e. for 
15 -N <>u,v sN) which is deducted from E c (u,v) to avoid numerical instabilities. 

Thus it can be seen that the similarity measure is modified by dividing the 
value of k by the size of the block W c . This is necessary so that the optimised value 
of k calculated for one image sequence can be used at all scales and resolutions (i.e. 
regardless of the size of the block W c chosen) for that sequence. 
20 The values of the response R c calculated using this equation are then used to 

calculate expected values of the velocity (u cc , and the corresponding covariance 
matrices using equations (4), (5) and (8) above. However, in this embodiment the 
calculation of the velocities (u^ v a )is further modified by using only candidate 
velocities which have probabilities above a certain threshold a in the velocity 
25 estimate of equations (4) and (5) however all candidate velocities are used in the 

covariance calculation. Thus in this embodiment the velocity estimates are calculated 
as follows:- 
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(20) 



T 

V 

cc 



where, 




0 otherwise 



R c (u,v) if R c (u,v)> a 



The threshold a is defined as follows: 



5 



a = m- T(m- m) with T e [0,l] 



where in and in are the maximum and minimum of the probability mass function 



R c (u,v) respectively. 

Thus it can be seen that if 7 is set to 1, the threshold a becomes the minimum 
value of R c , meaning that all values of the candidate velocities are used in the 

10 calculation, and the calculation becomes equivalent to that in the Singh approach. If 
7=0, on the other hand, the threshold a becomes the maximum value of the response 
so that only the candidate velocity with maximum probability is taken as the 
estimated velocity. Thus the estimate would be totally biased towards the 
predominant mode. In fact the value of T 9 optimised in the same optimisation 

15 process as that used for k, to be explained below, in practice will be between zero and 



The estimates of velocity and the co variance matrices are used together with 
neighbourhood information in the iterative process described above to calculate the 
optimised values of velocity in accordance with equation (12) above. 



together in a 2D space. In step 40 the sequence of images is taken and in step 41 the 



one. 



20 



Figure 4 illustrates schematically how the values of k and Tare optimised 
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values of k and T are initialised. Then in step 42 the image velocity is estimated 
using the initial values ofk and T. These initial values may be chosen from 
experience based on the type of imaging equipment and the subject of the imaging 
sequence. The process is relatively robust to the choice of k and T, so, for example, 

5 initial values of T=0.5 and A=0.5 may be suitable for an ultrasound imaging 

sequence. Having calculated the image velocity in step 42 it is then possible in step 
43 to register all of the subsequent frames to the first frame. "Registering" frames is 
equivalent to superimposing the images one upon the other and adjusting their 
relative position to get the best match. In practice the process involves correcting the 

10 subsequent frames for motion using the calculated image velocity. Having registered 
the frames a registration error £ is calculated using an error function in step 44. As 
an example, the error function may be a sum of square differences in the intensities 
of the frames. If the image velocity estimation were perfect, there would be no 
difference in intensities (as the motion correction would be perfect) and thus the error 

15 function would be zero. In practice, of course, the error function is non-zero and so 
in step 45 the values of k and T are varied to minimise the error function This may 
< be achieved using a multidimensional minimisation algorithm such as the Powell . 

algorithm (see William H. Press et al., "Numerical recipes in C: The art of scientific 
computing", Cambridge University Press). The optimisation process may be 

20 continued until the change in the value of the error function £ is below a certain 
threshold. In one experiment to compensate a breast compression sequence for 
distortion the optimal values were found to be T=0.660 and A=0.237. Figure 6 shows 
the results of the experiment conducted on the ultrasound breast data. The error 
shown is the registration error %. Two observations can be made: 

25 1 . For T = 0, the velocity estimation is equivalent to taking the argument 

of the maximum of the probability. Hence, theoretically, the parameter k does not 
have any influence on the result. This can easily be observed in this experiment, and 
it corresponds to the maximum error. In this case, the image velocity is quantified by 
the pixel resolution of the image, and hence the error on the image velocity is of the 

30 order of the pixel resolution. Furthermore, this approach is not robust against noise. 
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This explains this high error on the velocity estimation. 

2. For T=l (as in Singh), the velocity estimation is equivalent to taking the 
mean of the probability. The results are better than for T=0, but do not correspond to 
the optimal value. This result can be explained by the fact that taking the mean of the 
5 probability as an estimates of the velocity is not very precise and may lead to biased 
estimation if the pdf is not mono-modal or non-well-peaked pdfs. Observe as well the 
expected functional dependence between the two parameters (T and A:). This last 
point indicates that the search for the optimal values of T and k should be done in the 
2D space. In this experiment the optimal values are T= 0.660 and k = 0.237. Thus 

1 0 showing a clear distinction from the Singh result. 

It should be noted that the improvements above may be used in a coarse-to- 
fine strategy, i.e. a multiresolution approach in which velocities are first estimated at 
a low resolution, then at a next finer resolution these estimates are used as a first 
guess in the estimation process and the estimates are refined. This means that instead 

15 of searching in the window around y) in the second frame, one can search around 
(art-He,/ , y +v cy/) where u at and v^, are velocity estimates propagated from the coarser 
level. This approach is computationally efficient. Further, the image velocity 
estimation may be concentrated on regions of the image in motion, rather than 
conducted over the whole image. 

20 Figure 5 illustrates the process flow for the image velocity estimation given 

values of k and T(e.g. initial values if this is the first estimate for a given application, 
or optimised values). Firstly, in step 50, a sequence of image frames is taken. Then, 
in step 51, the similarity measure across three frame sets of the sequence is calculated 
using the CD 2 . bis similarity measure, i.e. using equation (18) at the desired scale and 

25 resolution. "Resolution" means whether one is sampling every pixel, or only certain 
pixels in the block W c and "scale" refers to how far the block is displaced in the 
search window W s9 e.g. by one pixel, or by several pixels. Having calculated the 
similarity values, the value of the response R c can be calculated in step 52 using 
equation (19). Then in step 53 the value of is calculated using equation (20) and 

30 the corresponding covariance matrix using equation (8). In step 54 the value of 
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U and the covariance for the neighbourhood estimate is calculated using equations 
(6), (7) and (9). Then in step 55 the conservation and neighbourhood information are 
fused using the iterative process of equation (12) to give an optimised velocity 
estimate U op . 

5 As indicated by step 56, the process may be repeated at finer scales and 

resolutions, with the computational burden being eased by making use of the image 
velocity estimate already obtained. 

The above improvements in the block matching technique are particularly 
successful in allowing tracking of cardiac boundary pixels in echocardiographic 
1 0 sequences. The block matching steps may be concentrated in a ribbon (band) around 
a contour defining the cardiac border to reduce the computational burden. However, 
the technique is applicable to other non-cardiac applications of ultrasound imaging. 
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CLAIMS 

1. A method of processing a sequence of image frames to estimate image 
velocity through the sequence comprising: 
5 block matching using a similarity measure by comparing the intensities in 

image blocks in two frames of the sequence and calculating the similarity between 
the said blocks on the basis of their intensities, calculating from the similarity a 
probability measure that the two compared blocks are the same, and estimating the 
image velocity based on the probability measure, wherein the probability measure is 
1 0 calculated using a parametric function of the similarity which is independent of 
position in the image frames. 

2. A method according to claim 1 wherein the parameters of the parametric function 
are independent of position in the image frames. 

15 

3. A method according to claim 2 wherein at least one of the parameters is optimised 
by coregistering the frames in the sequence on the basis of the calculated image 
velocity, calculating a registration error and varying at least one of the parameters to 
minimise the registration error. 

20 

4. A method according to claim 3 wherein the registration error is calculated from 
the differences of the intensities in the coregistered frames. 

5. A method according to claim 4 wherein the registration error is calculated from 
25 the sum of the squares of the differences of the intensities in the coregistered frames. 

6. A method according to any one of the preceding claims further comprising the 
step of normalising the calculated similarity with respect to the size of the block and 
calculating the probability measure on the basis of the normalised similarity. 
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7. A method according to claim 6 wherein the calculated similarity is normalised by 
dividing it by the number of image samples in the block. 

8. A method according to claim 6 wherein the calculated similarity is normalised by 
5 dividing it by the number of pixels in the block. 

9. A method according to any one of the preceding claims wherein the probability 
measure is a monotonic function of the similarity. 

10 10. A method according to any one of the preceding claims wherein the probability 
measure is thresholded such that motions in the image velocity whose probabilities 
have a predefined relationship with a threshold are ignored. 

1 1 . A method according to claim 10 wherein the threshold is optimised by 

15 coregistering the frames in the sequence on the basis of the calculated image velocity, 
calculating a registration error and varying the threshold to minimise the registration 
error. 

12. A method according to claim 10 or 1 1 wherein the threshold is positionally 
20 independent. 

13. A method according to claim 10, 1 1 or 12 wherein the threshold and parameters 
are optimised together. 

25 14. A method according to any one of the preceding claims further comprising 
normalising the intensities in the two blocks to have the same mean and standard 
deviation before calculating said similarity. 

15. A method according to any one of the preceding claims wherein the similarity 
30 measure is the CD 2Jbl5 similarity measure. 
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16. A method according to any one of the preceding claims wherein the block 
matching is conducted across three frames of the sequence by comparing the 
intensities in blocks in the first and third and the second and third of the three frames 

5 and calculating the similarity from said compared intensities. 

17. A method according to claim 16 wherein the blocks in the first and second 
frames are blocks calculated as corresponding to each other on the basis of a previous 
image velocity estimate. 

10 

18. A method of processing a sequence of image frames to estimate image velocity 
through the sequence comprising: 

block matching using a similarity measure by comparing the intensities in 
image blocks in three frames of the sequence by comparing the intensities in blocks 
1 5 in the first and third and the second and third of the three frames, and calculating the 
similarity between the said blocks on the basis of their intensities. 

19. A method according to claiml8 wherein the blocks in the first and second frames 
are blocks calculated as corresponding to each other on the basis of a previous image 

20 velocity estimate. 

20. A method according to claim 19 comprising defining for each block in the second 
frame a search window encompassing several blocks in the third frame, and 
calculating the similarity of each block in the search window to the said block in the 

25 second frame and to the corresponding position of the said block in the first frame 
based on the previous image velocity estimate. 

21. A method of processing a sequence of image frames to estimate image velocity 
through the sequence comprising: 

30 block matching using a similarity measure by comparing the intensities in 



.1 
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image blocks in two frames of the sequence and calculating the similarity between 
the said blocks on the basis of their intensities, further comprising normalising the 
intensities in the two blocks to have the same mean and standard deviation before 
calculating said similarity. 

5 

22. A method according to claim 21 wherein the similarity measure is the CD 2Jbis 
similarity measure. 

23. A method according to claim 21 or 22 wherein the block matching is conducted 
10 across three frames of the sequence by comparing the intensities in blocks in the first 

and third and the second and third of the three frames and calculating the similarity 
from said compared intensities. 

24. A method according to claim 23 wherein the blocks in the first and second 

1 5 frames are blocks calculated as corresponding to each other on the basis of a previous 
image velocity estimate. 

25. A method according to any one of the preceding claims wherein the image 
velocity estimate is refined by modifying the image velocity estimate at each position 

20 in the image with the estimated image velocity at surrounding positions. 

26. A method according to any one of the preceding claims wherein the images are 
medical images. 

25 27. A method according to any one of the preceding claims wherein the images are 
ultrasound images. 

28. Image processing apparatus comprising an image velocity estimator adapted to 
estimate image velocity in accordance with the method of any one of the preceding 
30 claims. 
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29. A computer program comprising program code means for executing on a 
programmed computer the method of any one of claims 1 to 27. 

30. A computer-readable storage medium storing a computer program according to 
5 claim 29. 
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Fig.4. 



Sequence of image frames 



Initialise /cand 7 
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Optical Flow estimation using 
/cand 7 
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Register all frames to the first 
frame 
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Calculate registration error 



45 
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Vary /cand 7 and repeat 
steps 42 to 44 to minimise 
registration error and obtain 
optimised values of /cand 7 
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Fig. 5. 



Input frame sequence 
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Calculate similarity across 
three frame sets of the 
sequence using the CD2_bis 
Eq 18 at the desired scale 
and resolution 
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Calculate Rc using Eq 19 
and optimised values of k 
and T for the application 
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Calculate lice using Eq 20 
and Sec using Eq 8 



54^ 



Initialise the neighbourhood 
information Ubar using Eq 6 
and 7 and Sn using Eq 9 
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Fuse the neighbourhood 
information using the iterative 
process of Eq. 1 2 to give 

Uop (update the 
neighbourhood information 
Ubar and Sn using the actual 
velocity estimate at every 
iteration 



56- 



Repeat at finer scales and 
resolutions if desired 
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